library(dplyr)
library(ggplot2)
library(readr)
library(viridis)
library(cowplot)
library(extrafont)
#
fte_theme <- function() {
    ## inspiré de http://minimaxir.com/2015/02/ggplot-tutorial/
    ## Generate the colors for the chart procedurally with RColorBrewer
    palette <- RColorBrewer::brewer.pal("Greys", n=9)
    color.background = palette[2]
    color.grid.major = palette[3]
    color.axis.text = palette[6]
    color.axis.title = palette[7]
    color.title = palette[9]
    ## Begin construction of chart
    theme_bw(base_size=9) +
        ## Set the entire chart region to a light gray color
        theme(
            text = element_text(family = "Ubuntu Light"),
            panel.background=element_rect(
                fill=alpha(color.background, 1), color=color.background),
            plot.background=element_rect(
                fill=alpha(color.background, 0.5), color=color.background),
            panel.border=element_rect(color=color.background),
            panel.grid.major=element_line(
                color=color.grid.major,size=.25, linetype = "dotted"),
            panel.grid.minor=element_blank(),
            axis.ticks=element_blank(),
            strip.text.y = element_text(angle = 0, size = 8),
            legend.position="none", # Format the legend, but hide by default
            legend.background = element_blank(),
            legend.text = element_text(size=9,color=color.axis.title),
            legend.key = element_rect(
                fill = scales::alpha("gray", 0),
                colour = scales::alpha("gray", 0)),
            ## Set title and axis labels, and format these and tick marks
            plot.title = element_text(
                color = color.title, size = 14, face = "bold", vjust = 1.25,
                hjust = 0, family = "Ubuntu"),
            axis.text.x = element_text(size = 7,color = color.axis.text),
            axis.text.y = element_text(size = 7,color = color.axis.text),
            axis.title.x = element_text(
                size = 8, color = color.axis.title, vjust = 0, hjust = 0.9),
            axis.title.y = element_text(
                size = 8, color = color.axis.title, vjust = 0.9, angle = 0 ),
            ## Plot margins
            plot.margin = unit(c(0.35, 0.2, 0.3, 0.35), "cm")
        )
}
#
theme_set(theme_bw() + fte_theme())
#
legend_position <- function(x=NULL, y=NULL) {
    custom_legend <- theme(
        legend.margin = unit(-0.3,"lines"),
        legend.key = element_rect(fill = scales::alpha("gray", 0),
                                  colour = scales::alpha("gray", 0)))
    if (is.null(x) & is.null(y)) {
        custom_legend + theme(legend.position = "bottom")
    } else {
        custom_legend + theme(legend.position = c(x, y))
    }
}
#
data_location <- "../../data/phruscle_snpcall.csv"
#
snp <- read_csv(data_location, col_types = "ccciciccid", trim_ws = TRUE) %>%
    mutate(
        name = gsub("-1073.+$", "", name),
        base = toupper(base)#,
    )
#
find_mutant <- function(name) {
    if      (grepl("ws", name)) "ws"
    else if (grepl("sw", name)) "sw"
    else if (grepl("W", name )) "w"
    else if (grepl("S", name )) "s"
    else ""
}
#
# neat little trick to reduce time of rowwise application of find_mutant.
snp <- snp %>%
    group_by(name) %>%
    summarise(count = n()) %>%
    rowwise() %>%
    mutate(mutant=find_mutant(name)) %>%
    inner_join(snp, .) %>%
    filter(mutant != "")

head(snp)

snp %>%
    ggplot(data = ., aes(x = qual )) +
    geom_histogram(binwidth = 1)

snp %>%
    ggplot(data = ., aes(x = refp, y = qual, color = qual )) +
    geom_point(alpha = 0.1, size = 0.1) +
    geom_line(aes(group = name), alpha = 1/20, size = 0.1 ) +
    scale_color_viridis()

snp %>%
    filter(qual < 30) %>%
    ggplot(data = ., aes(x = refp, y = qual )) +
    geom_point(alpha = 0.1, size = 0.1)

snp %>%
    group_by(name) %>%
    summarise(median = median(qual)) %>%
    qplot(data = . , median, binwidth = 0.1) +
    geom_vline(xintercept = 40, color = "red")

snp %>%
    group_by(name) %>%
    summarise(mean = mean(qual)) %>%
    qplot(data = . , mean, binwidth = 0.1) +
    geom_vline(xintercept = 40, color = "red")

keep_clean_only <- function(data, qual_ = 40) {

    get_dirty_seq <- function(data)
        group_by(data, name) %>%
            summarise(median = median(qual)) %>%
            filter(median < qual_) %>%
            select(name) %>% unlist()

    data %>%
        filter(!(name %in% get_dirty_seq(data)))
}

snp %>%
    group_by(name, mutant) %>%
    summarise(len = max(refp) - min(refp)) %>%
    ggplot(data = ., aes(x = len )) +
    geom_histogram(binwidth = 1) +
    facet_grid(mutant ~ .)

snp %>%
    filter(cons == "x" | cons == "X") %>%
    ## filter(refp < 750 & refp > 40) %>%
    ggplot(aes(x = refp)) +
    geom_histogram(binwidth = 1)

find_polarity <- function(reference, experimental) {
    is_w <- function(base) { ifelse(base == "A" | base == "T", TRUE, FALSE)}
    is_s <- function(base) { ifelse(base == "C" | base == "G", TRUE, FALSE)}

    if      (is_w(reference) & is_w(experimental)) "ww"
    else if (is_w(reference) & is_s(experimental)) "ws"
    else if (is_s(reference) & is_w(experimental)) "sw"
    else if (is_s(reference) & is_s(experimental)) "ss"
    else ""
}

plot_align <- function(data, mutant_) {

    sort_by_tract_length <- function(mut) {
        data %>%
            filter(mutant %in% mut) %>%
            group_by(name, mutant) %>%
            filter(cons == "x") %>%           # dernière position correspondant au donneur
            summarise(len = min(refp)) %>%
            ungroup() %>%
            arrange(len) %>%
            {.$name}
    }

    is_inside_conv <- function(data, mut) {
        data %>%
            filter(mutant %in% mut, cons == "x") %>%
            group_by(name) %>%
            summarise(bascul = min(refp)) %>%
            inner_join(data, ., by = "name") %>%
            keep_clean_only() %>%
            ## garde les positions après le point de switch qui correspondent au génotype parental.
            filter(refp > bascul, cons == "X" | cons == "x") %>%
            mutate(tract = TRUE) %>%
            select(name, bascul, refp, tract) %>%
            left_join(data, .)
    }

    is_restoration <- function(data, mut) {
        ## détermine si le SNP est une restauration du génotype parental. Autrement
        ## dit les SNP qui sont dans la trace de conversion mais qui sont X, *ie*
        ## des SNPs parentaux.
        data %>%
            filter(mutant %in% mut, cons == "X", tract == TRUE)
    }

    ## default color
    red    <- RColorBrewer::brewer.pal(n = 4, "Set1")[1]
    blue   <- RColorBrewer::brewer.pal(n = 4, "Set1")[2]
    green  <- RColorBrewer::brewer.pal(n = 4, "Set1")[3]
    violet <- RColorBrewer::brewer.pal(n = 4, "Set1")[4]

    data <-
        is_inside_conv(data, mutant_) %>%
        group_by(name, mutant) %>%
        filter(mutant %in% mutant_, cons == "x" | cons == "X" ) %>%
        rowwise() %>%
        mutate(sens = find_polarity(refb, seqb)) %>%
        ungroup() %>%
        keep_clean_only()

    data$tract[is.na(data$tract)] <- FALSE

    ggplot(data = data,
           aes(x = factor(refp),
               y = factor(name, levels = sort_by_tract_length(mut=mutant_)))) +
        geom_point(aes(color = tract, size = qual, alpha=qual)) +
        ## représente les cas complexes
        geom_point(data = is_restoration(data, mutant_),
                   aes(alpha = qual), color = green, size = 5) +
        ## représente la séquence du donneur
        geom_text(aes(label = snpb, x = factor(refp), y = -3), color = blue,
                  vjust = -0.5, size = 2, family = "Ubuntu Light") +
        annotate("text", x = 1.2, y = -0.8, label = "Donneur : ", color = blue, size = 2 ) +
        ## représente la séquence du receveur.
        geom_text(aes(label = refb, x = factor(refp),
                      y = length(sort_by_tract_length(mutant_)) + 10),
                  color = red, vjust = 9, size = 2, family = "Ubuntu Light") +
        annotate("text", x = 1.2, y = length(sort_by_tract_length(mutant_)) + 5,
                 label = "Receveur : ", color = red, size = 2 ) +
        scale_color_brewer(palette = "Set1", guide = "none") +
        scale_alpha( range=c(1/5, 0.8), guide=FALSE ) +
        scale_size(range = c(0.5, 2), breaks = c(10, 50),
                   labels = c("Faible", "Forte")) +
        labs(x = "", y = "", size = "Qualité",
             title = paste("Alignement pour la manip", toupper(mutant_))) +
        theme(legend.direction = "horizontal",
              ## legend.position  = c(0.8, -0.2),
              legend.margin = unit(0,"lines"),
              ## legend.key = element_rect(fill = scales::alpha("gray", 0),
              ##                           colour = scales::alpha("gray", 0)),
              panel.grid.major.y = element_line(size = 0.1, linetype = "dotted")) +
        legend_position(0.8, 0.98)
}


plot_qual <- function(data, mutant_) {
    red    <- RColorBrewer::brewer.pal(n = 4, "Set1")[1]
    blue   <- RColorBrewer::brewer.pal(n = 4, "Set1")[2]
    data %>%
        filter(mutant == mutant_) %>%
        ggplot(aes(x = refp, y = qual)) +
        geom_point(aes(color = qual), size = 1/10, alpha = 1/20) +
        geom_smooth(se = FALSE, color = red, linetype = "dotted") +
        scale_color_continuous(high = blue, low = red) +
        labs(x = "Position sur la référence", y = "")
}

plot_align_qual <- function(data, mutant_) {
    plot_grid(
        data %>% plot_align(mutant_),
        data %>% plot_qual(mutant_),
        ncol = 1, rel_heights = c(8.5, 1.5),
        align = 'v', labels = c("  1", "  2")
    )
}

snp %>% plot_align_qual("ws")

snp %>% plot_align_qual("sw")

snp %>% plot_align_qual("s")

snp %>% plot_align_qual("w")

plot_grid(
    snp %>% plot_align_qual("s"),
    snp %>% plot_align_qual("w"),
    snp %>% plot_align_qual("ws"),
    snp %>% plot_align_qual("sw"),
    labels = c("A", "B", "C", "D"), ncol = 2
)

## interactive only
## save_as_a3("test.pdf")

snp %>%
    filter(cons == "x") %>%
    group_by(name, mutant) %>%
    summarise(len = max(refp) - min(refp) ) %>%
    ggplot(aes(x = factor(mutant, levels = c("s", "w", "sw", "ws")),
               y = len)) +
    geom_point(alpha = 0.2) +
    stat_summary(fun.y = "mean", geom = "point", color = "red", shape = "|",
                 size = 7) +
    labs(x = "mutant", y = "Longueur\ndu tract") +
    coord_flip()

snp %>%
    filter(cons == "x") %>%
    group_by(name, mutant) %>%
    summarise(len = max(refp) - min(refp) ) %>%
    group_by(mutant) %>%
    summarise(mean_len = mean(len))

snp %>%
    filter(cons == "x") %>%
    group_by(name, mutant) %>%
    summarise(len = max(refp) - min(refp) ) %>%
    ggplot(data = ., aes(x = len, fill = mutant )) +
    geom_histogram(binwidth = 5, position = "dodge") +
    facet_grid(mutant ~ .) +
    scale_fill_brewer(palette = "Set1") +
    labs(x = "Longueur de trace de conversion",
         y = "",
         title = "Distribution des\nlongueurs de trace de conversion")

get_gc_content <- function(data, mut=NULL, tidy = TRUE) {
    clean_seq <- function(x) toString(x) %>% gsub(", ", "", . )

    get_gc_rate <- function(dna_seq) {

    }
    if (is.null(mut)) mut <- c("ws", "sw", "w", "s")

    gc_content <- data %>% filter(mutant %in% mut) %>% group_by(name, mutant) %>%
        summarise(exp = seqinr::GC(seqb),
                  snp = seqinr::GC(snpb),
                  ref = seqinr::GC(refb))
    if (tidy) gc_content %>% tidyr::gather("type", "GC", 3:5)
    else gc_content
}


snp %>%
    get_gc_content() %>%
    filter(name != "psw21") %>%
    mutate(mutant = factor(mutant, levels = c("s", "w", "sw", "ws"))) %>%
    ## group_by(mutant, type) %>%
    ## summarise(gc_mean = mean(GC)) %>%
    ggplot(aes(y = GC,
               x = factor(type, levels = c("snp", "exp", "ref"),
                          labels = c("Donneur", "Expérimental", "Receveur")))) +
    geom_jitter(aes(color = mutant), width = 1/4, alpha = 1/3, size = 1) +
    geom_line(aes(group = name, color = mutant), alpha = 1/20) +
    stat_summary(fun.y = mean, geom = "point", color = "black", shape = "¦",
                 size = 6) +
    facet_grid(mutant ~ .) +
    labs(y = "% GC", x = "", title = "Comparaison des taux de GC") +
    scale_color_brewer(palette = "Set1") +
    coord_flip()

is_inside_conv <- function(data, mut) {
  data %>%
      filter(mutant %in% mut, cons == "x") %>%
      group_by(name) %>%
      summarise(bascul = min(refp)) %>%
      inner_join(data, ., by = "name") %>%
      keep_clean_only() %>%
      ## garde les positions après le point de switch qui correspondent au génotype parental.
      filter(refp > bascul, cons == "X" | cons == "x") %>%
      mutate(tract = TRUE) %>%
      select(name, bascul, refp, tract) %>%
      left_join(data, .)
}

## wip
snp %>%
    is_inside_conv(c("ws", "sw", "s", "w")) %>%
    filter(name == "pS10")

## /wip



sessionInfo()